In this project, I will be investigating and forecasting electricity consumption in the UK. I will be comparing two different, both very popular forecasting methodologies: Holt-Winters exponential smoothing and ARIMA models. I will be comparing and evaluating both methodologies against the naive model which simply forecasts each observation to be equal to the previous one (Y_t = Y_t-1). I will be using four different performance metrics:
1.) Mean Squared error (MSE): MSE simply takes the difference between the forecasted value Yhat and the actual value Y, squares it for each observation and then takes the average over all observations. The reason for this is so that negative and positive values don’t cancel each other out. Additionally, this also penalises predictions that are very far off more than ones that are only slighty wrong. 2.) Mean Absolute error (MAE): MAE also computes the difference between forecasted value and actual value, but instead of squaring it simply takes the absolute value of the difference. It then computes the average among all the differences. This method is preferable if we want all errors to be weighted equally, and don’t want errors that are further off to be penalised higher.
3.) Root Mean Squared Error (RMSE): RMSE is simply the square root of the MSE. The reason for this is so the error is on the same scale as the original data, giving better insight into how the error term compares to the data.
4.) Mean Absolute Percentage Error (MAPE): MAPE is the average ratio between the absolute error and the absolute value of the actual value (Y).
Let’s take a look at the data:
data(USgas)
ts_info(USgas)
## The USgas series is a ts object with 1 variable and 238 observations
## Frequency: 12
## Start time: 2000 1
## End time: 2019 10
ts_plot(USgas, Xtitle = "Year",
Ytitle = "Natural gas consumption in Million cubic feet",
title = "Monthly US Natural gas consumption")
Looks like there is quite a strong seasonal component in the data, which makes sense as a lot of heating is powered using natural gas so there would naturally be more demand for gas during the winter months compared to the summer months. We can also see a slightly increasing trend from about 2010 onwards. Let’s take a look at a decomposition of the time series and see if this confirms these first visual impressions:
decompose(USgas) %>% plot()
Decomposition reveals a trend as well, beginning around 2005 but becoming more pronounced from about 2010 onwards. I’ll be splitting the data into training and test set, the test set will be the last 12 months of the data. I’ll start off by fitting a naive model, which simply uses the average of all preceding observations. This will be a good benchmark to compare the performance of more complex models to.
USgas_par <- ts_split(USgas, 12)
train <- USgas_par$train
test <- USgas_par$test
naive_mod <- naive(train,h=12)
fc_naive <- forecast(naive_mod,h=12)
accuracy(fc_naive, test)
## ME RMSE MAE MPE MAPE MASE
## Training set -1.028444 285.6607 228.5084 -0.9218463 10.97123 1.984522
## Test set 301.891667 499.6914 379.1417 9.6798015 13.28187 3.292723
## ACF1 Theil's U
## Training set 0.3761105 NA
## Test set 0.7002486 1.499679
test_forecast(actual = USgas,
forecast.obj = fc_naive,
test = test)
As we can see, the naive model simply takes the average of all previous observations without any considerations of seasonal or trend data. Let’s fit a Holt Winters exponential smoothing model instead and see how well this captures the data and whether we get an improvement over the naive model:
hw_mod <- HoltWinters(train)
hw_mod
## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = train)
##
## Smoothing parameters:
## alpha: 0.371213
## beta : 0
## gamma: 0.4422456
##
## Coefficients:
## [,1]
## a 2491.9930104
## b -0.1287005
## s1 32.0972651
## s2 597.1088003
## s3 834.9628439
## s4 353.8593860
## s5 318.1927058
## s6 -173.0721496
## s7 -346.6229990
## s8 -329.7169608
## s9 -112.1664217
## s10 -140.3186476
## s11 -324.5343787
## s12 -243.9334551
We can see that the model takes into account the average level of the preceding observations (alpha = 0.37) as well as a strong seasonal component (gamma = 0.44). It does not, however, take into account a trend (beta = 0).
fc_hw <- forecast(hw_mod, h = 12)
accuracy(fc_hw, test)
## ME RMSE MAE MPE MAPE MASE
## Training set 7.828361 115.2062 87.67420 0.2991952 4.248131 0.7614222
## Test set 51.013877 115.1555 98.06531 1.7994297 3.766099 0.8516656
## ACF1 Theil's U
## Training set 0.21911103 NA
## Test set -0.01991923 0.3652142
This is a significant improvement over the naive model, with RMSE reducing by more than 75% from 500 to only 115. We can see a similar reduction when looking at MAPE and MAE.
test_forecast(actual = USgas,
forecast.obj = fc_hw,
test = test)
Looking visually at our predictions, this confirms what the error metrics tell us. The model beautifully captures the seasonal patterns.
Let’s see how a different class of models, the ARIMA model, performs for this dataset. Arima models rely on a combination of two modelling techniques, the AR and MA processes. The AR process basically explains a model’s future values as a linear combination of it’s previous values, or lags. It requires the time series to be stationary, meaning it can have an increasing or decreasing trend or variance. Lots of time series are not stationary however, be it due to trends, varying seasonal patterns or a certain number of random events. Fortunately, there is a way to deal with these time series as well using a technique called differencing. This mean subtracting the value from one cycle prior (so for example in a daily series, Yt = Yt - Yt-365). After doing this, we simply analyse the difference between these two time points. Forecasts made on these differences easily can be converted back to an actual time series.
Due to the strong seasonality in the data, a SARIMA model, which has three additional parameters to capture seasonality, will likely be the most appropriate. To get a better understanding of the serial correlation of the data (how much each observation is correlated with previous observations, or lags), I will plot the ACF and PACF plots.
par(mfrow=c(1,2))
acf(train, lag.max = 60)
pacf(train, lag.max = 60)
The ACF plot indicates that there is a strong positive correlation with previous observations of the same season, and a weaker but still significant negative correlation with observations of the opposite season. All of this does not come as a surprise since we know that natural gas consumption is highly seasonal. We can also see that the correlation is decaying over time, indicating that the series is in fact not stationary and we will need to do some differencing to continue with our ARIMA model.
USgas_d12 <- diff(train, 12) #specifying the period of differencing, in this case 12 months
ts_plot(USgas_d12,
title = "US Monthly Natural Gas consumption - First Seasonal
Difference",
Ytitle = "Billion Cubic Feet (First Difference)",
Xtitle = "Year")
This yields the following. We have successfully removed the trend in the series, however there still exists non-constant variance which we will need to account for by differencing once more:
USgas_d12_1 <- diff(diff(USgas_d12, 1))
ts_plot(USgas_d12_1,
title = "US Monthly Natural Gas consumption - First Seasonal and
Non-Seasonal Differencing",
Ytitle = "Billion Cubic Feet (Difference)",
Xtitle = "Year")
This looks better, let’s look again at the ACF and PACF plots after the transformation:
par(mfrow=c(1,2))
acf(USgas_d12_1, lag.max = 60)
pacf(USgas_d12_1, lag.max = 60)
As we can see in the plot, after the transformations the correlation with the lags seem to be tailing off very quickly. Let’s try fit an arima model using auto.arima. This automatically determines the best parameters for the AR and MA processes, as well as doing the required differencing for us.
USgas_arima_mod <- auto.arima(train)
USgas_arima_mod
## Series: train
## ARIMA(2,1,1)(2,1,1)[12]
##
## Coefficients:
## ar1 ar2 ma1 sar1 sar2 sma1
## 0.4301 -0.0372 -0.9098 0.0118 -0.2673 -0.7431
## s.e. 0.0795 0.0755 0.0452 0.0886 0.0826 0.0748
##
## sigma^2 estimated as 10446: log likelihood=-1292.83
## AIC=2599.67 AICc=2600.22 BIC=2623.2
USgas_test_fc <- forecast(USgas_arima_mod, h = 12)
accuracy(USgas_test_fc, test)
## ME RMSE MAE MPE MAPE MASE
## Training set 5.846395 97.81607 73.42635 0.1171733 3.522336 0.6376842
## Test set 37.858513 103.23333 81.47035 1.3112392 3.261787 0.7075437
## ACF1 Theil's U
## Training set -0.004151171 NA
## Test set -0.046726335 0.3404228
Not bad, we reduced RMSE from 115 using HoltWinters to just 103. Let’s take a look at the plot of the forecast:
test_forecast(USgas,
forecast.obj = USgas_test_fc,
test = test)
This model also captures the seasonal component very well, while also including the trend. We can see that the peak for the forecasted time points is just a little bit higher than in the previous years, which is what we would expect looking at previous year’s changes. This is something the Holt Winters model did not pick up on, which is reflected in the slightly lower RMSE. To check how well the model fits our data, let’s perform a residuals check:
checkresiduals(USgas_arima_mod)
##
## Ljung-Box test
##
## data: Residuals from ARIMA(2,1,1)(2,1,1)[12]
## Q* = 24.949, df = 18, p-value = 0.1263
##
## Model df: 6. Total lags used: 24
Everything looks good here, we can see that the residuals resemble white noise and are evenly distributed with a mean of zero. There are no significant autocorrelations in the lags, and the Ljung-Box test also does not reject the null-hypothesis of no auto correlation with a non-significant p-value of 0.13.
And there we go, we’ve successfully fitted a few different models, both of which do a good job of capturing the data, but I think it’s fair to say that the ARIMA model has come out as the clear winner in this case. Let’s take another look at it’s predictions for the last 12 months, complete with a confidence interval:
plot_forecast(USgas_test_fc,
title = "US Natural Gas Consumption - Forecast",
Ytitle = "Billion Cubic Feet",
Xtitle = "Year")